/* Multivariate Poisson-lognormal VAR(p) model */

# Model definition
model { for (i in 4:N) 
      	{ 
	# Define the random effects
	eta[i,1:m] ~ dmnorm(nought, P)

	# Define the Poisson RVs and the mean function
	SCcites[i] ~ dpois(mu[1,i])
	ACcites[i] ~ dpois(mu[2,i])
	DCcites[i] ~ dpois(mu[3,i])

	# Lag dynamics

	mu[1,i] <- SC1_SC*SCcites[i-1] + AC1_SC*ACcites[i-1] + DC1_SC*DCcites[i-1] +
		   exp(Intercepts[1] + age_SC*age[i] + distance_SC*distance[i] + overruled_SC*overruled[i] +
		   eta[i,1]);

	mu[2,i] <- SC1_AC*SCcites[i-1] + AC1_AC*ACcites[i-1] +  DC1_AC*DCcites[i-1] +
		   exp(Intercepts[2] + age_AC*age[i] + distance_AC*distance[i] + overruled_AC*overruled[i] +
		   eta[i,2]);

	mu[3,i] <- SC1_DC*SCcites[i-1] + AC1_DC*ACcites[i-1] +  DC1_DC*DCcites[i-1] +
		   exp(Intercepts[3] + age_DC*age[i] + distance_DC*distance[i] + overruled_DC*overruled[i] +
		   eta[i,3]);

	}

	# Cov mtx and correlations

    	Sigma <-inverse(P)
	for (i in 1:m) {
	    for (j in 1:m) {
	    C[i,j] <- Sigma[i,j]/sqrt(Sigma[i,i]*Sigma[j,j])}}

	# Priors: precision matrix, random effects, correlations, coefficients.

	P ~ dwish(R, m)	
	for (i in 1:m) 
	    { nought[i] <- 0
	    	for (j in 1:m) 
		{ R[i,j] <- equals(i,j)
	    	}
	    }

	# Prior for the intercepts -- SZ style from the moments
		
	Intercepts[1] ~ dnorm(0, 0.3)
	Intercepts[2] ~ dnorm(0, 0.3)
	Intercepts[3] ~ dnorm(0, 0.3)

 	# Prior for the AR coefficents -- SZ style from the moments

	# First lag

  	SC1_SC ~ dnorm(1, 1)
  	AC1_AC ~ dnorm(1, 1)
  	DC1_DC ~ dnorm(1, 1)

  	AC1_SC ~ dnorm(0, 2)
  	DC1_SC ~ dnorm(0, 2)
  	SC1_AC ~ dnorm(0, 2)
  	DC1_AC ~ dnorm(0, 2)
  	SC1_DC ~ dnorm(0, 2)
  	AC1_DC ~ dnorm(0, 2)

	# Prior for exogenous variables
	age_SC ~ dnorm(0, 0.6)
	age_AC ~ dnorm(0, 0.6)
	age_DC ~ dnorm(0, 0.6)
	distance_SC ~ dnorm(0, 0.6)
	distance_AC ~ dnorm(0, 0.6)
	distance_DC ~ dnorm(0, 0.6)
	overruled_SC ~ dnorm(0, 0.6)
	overruled_AC ~ dnorm(0, 0.6)
	overruled_DC ~ dnorm(0, 0.6)
	
}
